Beyond Quit[]: 
Master Your Memory in the Wolfram Language

Anthony Zupnik
anthonyz@wolfram.com

Technical Content Developer

Why listen to this talk?

Running out of memory

Slow running code

Slow running dynamics

Mathematica locking up when executing code

Incentives!

Computation History

Quit the Kernal for a clean slate

Quit[]

How much memory is in use

MemoryInUse[]

Find out start memory, add data and compute the difference.

(start = MemoryInUse[];
  data = RandomReal[1, {1000, 1000}];
  Print[MemoryInUse[] - start];)

Lets have a closer look at assigning values

MemoryInUse[]

m = RandomReal[1, {1000, 1000}];

MemoryInUse[]

ByteCount[m]

Clear[m]; MemoryInUse[]

Downvalues[Out] will give all the Downvalues associated with Out.

DownValues[Out][[1 ;; 4]]

ByteCount /@ DownValues[Out]

In Downvalues can be seen in the same way

ByteCount /@ DownValues[In]

Run this a few times

ByteCount /@ DownValues[Out]

MemoryInUse[]

To Remove all Output you can do:

Unprotect[Out]; Clear[Out]

MemoryInUse[]

To Remove all In you can do:

Unprotect[In]; Clear[In]

MemoryInUse[]

You can however remove a specific Output (Line is not reset so we use $Line)

$Line

$Line = 1

RandomReal[1, {1, 1}]

Unprotect[Out]; Out[2] =.

Out[2]

m = RandomReal[1, {1000, 1000}];

ByteCount[m]

Clear[m]; MemoryInUse[]

Unprotect[In]; Clear[In]

Unprotect[In]; Clear[Out]

MemoryInUse[]

m = RandomReal[1, {1000, 1000}];

MemoryInUse[]

Conviniently there is also a function to remove all Ins and Outs

<<Utilities`CleanSlate`

Utilities`CleanSlate`ClearInOut[]

DownValues[Utilities`CleanSlate`ClearInOut]

History Length

Quit[]

$HistoryLength

MemoryInUse[]

RandomReal[1, {1000, 1000}];

MemoryInUse[]

$HistoryLength = 0

RandomReal[1, {1000, 1000}];

MemoryInUse[]

Quit[]

MemoryInUse[]

$HistoryLength = 0

RandomReal[1, {1000, 1000}];

MemoryInUse[]

{DownValues[In], DownValues[Out]}

You can store something like this in init.m if you want it to initialise each and every time you run Mathematica

$UserBaseDirectory // SystemOpen

Quit[]

$HistoryLength

Special Data Types

Mathematica will often optimise your data structures from lists of lists into more memory efficient structures.

PackedArray

PackedArray is useful for where all elements are of the same type (Integer or Real), it is used automatically behind the scenes - and may be checked using Developer`PackedArrayQ

realVec = Table[RandomReal[], {100000}];
intVec = Table[RandomInteger[], {100000}];

{Developer`PackedArrayQ[realVec], Developer`PackedArrayQ[realVec]}

Note that as the type of each element does not need to be tracked it is a much smaller object:

ByteCount[realVec]

The structure is destroyed when an element is added with Append as Mathematica does not check if all elements are of the same type.

ByteCount[Append[realVec, 1]]

Developer`ToPackedArray[expr, type] constructs a PackedArray of type from the evaluation of expr.

appendedRealVec = Developer`ToPackedArray[Append[realVec, 1], Real];
ByteCount[appendedRealVec]

Note that these PackedArrays are more efficient to compute with as a result of their compressed structure:

Total[realVec]; // AbsoluteTiming

Total[Developer`FromPackedArray[realVec]] // AbsoluteTiming

SparseArray

SparseArray is an efficient object for storing data where most elements are identical, typically they are considered zero, and only non-identical elements (and their positions) are specified.

SparseArray[{position -> value, ...}, {d1, d2, ..., dn}, value] generates a SparseArray with dimensions d1×d2×...×dn with all elements except those at position set as value.

sparseA = SparseArray[{{1, 1} -> a, {3, 8} -> 3, {8, 3} -> c}, {10, 10}, "empty"];
MatrixForm[sparseA]

Normal can be used to convert a SparseArray to an ordinary dense matrix, note again that the optimised data structure is significantly smaller.

{ByteCount[sparseA], ByteCount@Normal[sparseA]}

ByteCount@Normal[sparseA]

As an example of this in force, consider calculating the modulus of a 1000 × 10000 matrix with and without SparseArray

constArray = ConstantArray[0, {1000, 1000}];
Do[constArray[[RandomInteger[{1, 1000}], RandomInteger[{1, 1000}]]] = RandomReal[100], {400}];
ByteCount[constArray]

sparseArray = SparseArray[{}, {1000, 1000}];
Do[sparseArray[[RandomInteger[{1, 1000}], RandomInteger[{1, 1000}]]] = RandomReal[100], {400}];
ByteCount[sparseArray]

Sqrt[constArray.constArray]; // AbsoluteTiming

Sqrt[sparseArray.sparseArray]; // AbsoluteTiming

Plot Effeciency

MaxPlotPoints

Mesh to see plot points

{#, ByteCount[#]} &@ListLinePlot[Table[{x, Sin[x]}, {x, 0, 2 Pi, 0.1}], Mesh -> All, ImageSize -> 500]

{#, ByteCount[#]} &@ListLinePlot[Table[{x, Sin[x]}, {x, 0, 2 Pi, 0.1}], Mesh -> All, MaxPlotPoints -> 25, ImageSize -> 500]



Why is this important?

{#, ByteCount[#]} &@ListLinePlot[Sort[RandomVariate[NormalDistribution[0, 1], {1000000}]], ImageSize -> 300]

{#, ByteCount[#]} &@ListLinePlot[Sort[RandomVariate[NormalDistribution[0, 1], {1000000}]], MaxPlotPoints -> 500, ImageSize -> 300]

Consider a plot, allow Mathematica to automatical select a method for you

{ByteCount[#],#}&@Plot3D[(1/(1+Abs[-1+(x^2)+(y^2)])),{x,-8,8},{y,-8,8},PlotRange->All,ImageSize->500,Mesh->All]

MaxRecursion defines how many recursive subdivisions which can be made

{ByteCount[#],#}&@Plot3D[(1/(1+Abs[-1+(x^2)+(y^2)])),{x,-8,8},{y,-8,8},PlotRange->All,PlotPoints->150,MaxRecursion->0,Mesh->All,ImageSize->500]

{ByteCount[#],#}&@Plot3D[(1/(1+Abs[-1+(x^2)+(y^2)])),{x,-8,8},{y,-8,8},PlotRange->All,MaxRecursion->10,PlotPoints->10,Mesh->All,ImageSize->500]

Operating Dymanics

Automatically Mathematica will choose to make the performance based on speed when controller is active, and quality when released.

Manipulate[Plot3D[Im[Zeta[x + I y]], {x, -n, 0}, {y, -n, n}], {n, 1, 20}]

ControlActive changes which expresion used based on the the state of the controller, so we can control outselves what the quality is.

Manipulate[Plot3D[Im[Zeta[x + I y]], {x, -n, 0}, {y, -n, n}, PlotPoints -> ControlActive[10, 20], MaxRecursion -> ControlActive[1, 3]], {n, 1, 20}]

ContiniousAction determines if something is continuously evaluated or not

Manipulate[Plot3D[Im[Zeta[x + I y]], {x, -n, 0}, {y, -n, n}], {n, 1, 20, ContinuousAction -> False}]

Allow Mathematica to do the work so you don’t have to by using PerformanceGoal instead, although the below is done automatically

Manipulate[Plot3D[Im[Zeta[x + I y]], {x, -n, 0}, {y, -n, n}, PerformanceGoal -> ControlActive["Speed", "Quality"]], {n, 1, 20}]

If we care about quality then set to “Quality”

Manipulate[Plot3D[Im[Zeta[x + I y]], {x, -n, 0}, {y, -n, n}, PerformanceGoal -> "Quality"], {n, 1, 20}]

Additional Note on Dynamics

Dynamic Objects only evaluate when they are on screen.

x = 0;

Dynamic[x = x + 1]

Dynamic[x]

This is an expensive manipulate, it is best to hide it after use to preserve memory usages.

OpenerView[{"A Manipulate", Manipulate[
      (*If the number of points has been reduced, discard  points*)
      If[atomcount < Length[contents], contents = Take[contents, atomcount]];
      (*If the number of points has been increased, generate some random points*)
      If[atomcount > Length[contents], contents = Join[contents, Table[{RandomReal[{-size, size}, {2}], {Cos[#], Sin[#]} &[RandomReal[{0, 2 π}]]}, {atomcount - Length[contents]}]]];
      
      Grid[{{Graphics[{PointSize[0.02],
                (*Draw the container*)
                Line[size {{-1, -1}, {1, -1}, {1, 1}, {-1, 1}, {-1, -1}}],
                Blend[{Blue, Red}, temperature/0.3],
                Point[
                  
                  (*Start the main dynamic actions*)
                  Dynamic[
                    (*Reset the collision counter*)
                    collision = 0;
                    (*Check for mouse interaction and add points if there has been one*)
                    Refresh[If[pt =!= lastpt, atomcount++; AppendTo[contents, {pt, {Cos[#], Sin[#]} &[RandomReal[{0, 2 π}]]}]; lastpt = pt], TrackedSymbols -> {pt}];
                    (*Update the position of the points using their velocity values*)
                    contents = Map[{#[[1]] + #[[2]] temperature, #[[2]]} &, contents];
                    (*Check for and fix points that have exceeded the box in Y direction, incrementing the collision counter for each one*)
                    contents = Map[
                        If[Abs[#[[1, 2]]] > size, collision++; {{#[[1, 1]], 2 size Sign[#[[1, 2]]] - #[[1, 2]]}, {1, -1} #[[2]]}, #] &,
                        contents];
                    (*Check for and fix points that have exceeded the box in X direction, incrementing the collision counter for each one*)
                    contents = Map[
                        If[Abs[#[[1, 1]]] > size, collision++; {{2 size Sign[#[[1, 1]]] - #[[1, 1]], #[[1, 2]]}, {-1, 1} #[[2]]}, #] &,
                        contents];
                    pressurehistory = Take[PadLeft[Append[pressurehistory, collision/size], 200], 200];
                    Map[First, contents]]]},
              
              PlotRange -> {{-1.01, 1.01}, {-1.01, 1.01}}, ImageSize -> {250, 250}],
            
            (*Show the smoothed historical values*)
            Dynamic@Show[
                ListPlot[
                  Take[MovingAverage[pressurehistory, smooth], -100],
                  Joined -> True, ImageSize -> {250, 250}, AspectRatio -> 1, PlotLabel -> "pressure history", AxesLabel -> {"time", "pressure"}, PlotRange -> {0, Max[Max[pressurehistory], 1]}],
                If[showtheoritical,
                  Plot[0.634 atomcount temperature/size^2, {x, 0, 100}],
                  Graphics[]]
                ]}}],
      
      {{pt, {0, 1}}, {-1, -1}, {1, 1}, Locator, Appearance -> None},
      {{atomcount, 50, "number of atoms"}, 0, 300, 1},
      {{temperature, 0.05, "temperature"}, 0.002, 0.3},
      {{size, 1, "container size"}, 0.2, 1},
      {{smooth, 1, "measurement time"}, 1, 100, 1},
      {{showtheoritical, False, "show theoretical value"}, {True, False}},
      {pressurehistory, {{}}, ControlType -> None},
      
      {contents, {{}}, ControlType -> None},
      {lastpt, {{0, 1}}, ControlType -> None},
      SynchronousUpdating -> True, AutorunSequencing -> {2}, ControlPlacement -> Left]}]

Additional Memory Related Functions

Clearing the System Cache

Quit[]

MemoryInUse[]

Reevaluating may take a shorter amount of time based on what is cached by the system

Timing[N[Pi, 10^5];]

Timing[N[Pi, 10^5];]

MemoryInUse[]

ClearSystemCache[]

MemoryInUse[]

Timing[N[Pi, 10^5];]

Other Notes

Quit[]

We can find out how much CPU time is used in a session

TimeUsed[]

Maximum Memory Used at any given time in a session

MaxMemoryUsed[]

Run a command until the memory limit is exceeded

MemoryConstrained[Range[10^6], 10000]

Similar but on evaluation time

TimeConstrained[Do[i^2, {i, 1000000}], 0.5]

Lets take our list line plot from before, we could add some kind of constraint, either Memory or Time related.

Quit[]

Dynamic[Refresh[MemoryInUse[], UpdateInterval -> 0.1]]

start = MemoryInUse[]

If[(plot = MemoryConstrained[ListLinePlot[Sort[RandomVariate[NormalDistribution[0, 1], {1000000}]]], 100000000]) === $Aborted, plot = ListLinePlot[Sort[RandomVariate[NormalDistribution[0, 1], {1000000}]], MaxPlotPoints -> 500], plot]

MaxMemoryUsed[] - start

If[TimeConstrained[ListLinePlot[Sort[RandomVariate[NormalDistribution[0, 1], {1000000}]]], 2] === $Aborted, ListLinePlot[Sort[RandomVariate[NormalDistribution[0, 1], {1000000}]], MaxPlotPoints -> 500]]

Use DumpSave (Storing your data)

We can use dumpsave to save definitions of any symbols we would like. This is stored in binary.

Quit[]

x = 2;
y = 6;
z = 7;
a = Plot[Cos[x], {x, 1, 10}];

?Global`*

Plot[Sin[x], {x, 1, 10}]

SetDirectory[$TemporaryDirectory]

DumpSave["f.mx", {"Global`", Out, In}]

Quit[]

?Global`*

Note: All previously assigned definitions will be destroyed when you Get the .mx file

<< f.mx

We can see everything has been imported back.

?Global`*

DownValues[In]

DownValues[Out]

Importing/Exporting Different Types

SetDirectory[NotebookDirectory[]];

exportData=RandomInteger[{0,1000},{1000,100}];

Export["testData.xlsx",exportData];//AbsoluteTiming

Export["testData.csv",exportData];//AbsoluteTiming

DumpSave["testData.mx",exportData];//AbsoluteTiming

Quit[]

SetDirectory[NotebookDirectory[]];

Import["testData.xlsx"];//AbsoluteTiming

Import["testData.csv"];//AbsoluteTiming

Import[ "testData.mx"]//AbsoluteTiming

Expand memory usage

Launch Kernels

LaunchKernels[]

Mathematica supports both CPU and GPU-based parallelisation, using CUDA and OpenCL hardware when available.

Parallelize and its Relatives

Not all functions that can be parallelised know about it, there’s a neat function called Parallelize that will try to parallelise your code across all available CPU cores to you:

This searches for all Mersenne Primes over the range 1000 to 3000 but is relatively slow

Select[Range[1000, 3000], PrimeQ[2^# - 1] &] // AbsoluteTiming

Parallelize rapidly speeds this up. Note that the Workbench Profiler will give you information about what exactly is generated by Parallelize.

Parallelize[Select[Range[1000, 3000], PrimeQ[2^# - 1] &]] // AbsoluteTiming

In general Mathematica will always distribute definitions of functions (and variables) to all CPU kernels, for instance:

sqrFunc[x_] := x^2 + 1;
Parallelize[Map[(Pause[1]; sqrFunc[#]) &, Range[10]]] // AbsoluteTiming

Remember it isn’t always a good idea

data = RandomReal[1, {500, 100, 100}];

Map[Transpose, data]; // AbsoluteTiming

ParallelMap[Transpose, data]; // AbsoluteTiming

Reap and Sow

In general Mathematica functions operate in an immutable fashion, however there are many different strategies for appending new elements to a data structure.

AppendTo[sequence, elem] appends elem to sequence with the definition of sequence being updated accordingly.

randList = {1, 2, 3};
AppendTo[randList, RandomReal[]];
randList

Suppose that a great many elements need to be appended, according to some algorithm, it appears natural to use a Do loop:

SeedRandom[1];
AbsoluteTiming[appendToDo = Do[AppendTo[randList, RandomReal[x]], {x, 0, 40000}];]

This is extraudinarily slow, the value of appendToDo must be updated after each iteration making this both system and memory intensive.

One faster method would be to use Sow and Reap - Reap collects everything that is provided to Sow:

AbsoluteTiming[reapAndSow = Reap[Do[Sow[RandomReal[x]], {x, 0, 40000}]][[2, 1]];
  AppendTo[randList, reapAndSow];]

However, in this instance it is better to remember the second tip - use the minimal number of built in functions.

AbsoluteTiming[
  AppendTo[randList, Table[RandomReal[i], {i, 0, 40000}]];]

Process Data a bit at the time

Quit[]

SetDirectory[NotebookDirectory[]]

Open read allows us to read through line by line, word by word, etc to process data a bit at a time

readData = OpenRead["tb.csv"]

Read a Record through with Read[_,Record]

test = Read[readData, Record]

The line is a string, lets take the 5 and 6 parts.

StringSplit[test, ","][[{5,6}]]

We can write our processed results to a file if we wish, no argument will result in a temporary file being made.

writeData = OpenWrite["test.txt"]

Do[Write[writeData, StringSplit[Read[readData, Record], ","][[5 ;; 6]]], {100}]

Close[writeData]

Import["test.txt"]

openApp = OpenAppend["test.txt"]

Do[Write[openApp, StringSplit[Read[readData, Record], ","][[5 ;; 6]]], {100}]

Close[openApp]

% // SystemOpen

We can reset, or chose a position in the stream. Here we’re resetting to zero

SetStreamPosition[readData, 0]

Clear[x]

Now run a sow reap for more efficiency and read through all the files.

(data = Reap[While[(x=Read[readData, Record])=!= EndOfFile, Sow[ToExpression[StringSplit[x, ","][[5 ;; 6]]]]]][[2,1]]) // AbsoluteTiming

MaxMemoryUsed[]

Close[readData]

As a direct comparison to import

Import["TB.csv"][[All, {5, 6}]] // AbsoluteTiming

MaxMemoryUsed[]

Compile Low Level Code

The Wolfram Language is in general loosely typed, providing a rich programming environment with very flexiblity functionality. However, where efficiency of code is required compilation becomes important and so with it flexibility less so.

Mathematica allows most low-level functions to be compiled to either optimised Wolfram Language byte code or symbolically generared C Code.

Compiled Wolfram Language Code

The following low-level function is a prime candidate for compilation,

normalFunc = Function[{x},
      Block[{sum = 1.0, inc = 1.0}, Do[inc = inc*x/i; sum = sum + inc, {i, 10000}]; sum]];

Mapping it across Range[-50., 50., 0.25] is relatively slow:

arg = Range[-50., 50., 0.5];
AbsoluteTiming[Map[normalFunc, arg];]

Compile[{arg1, type}, {arg2, type}, ...}, expr] creates a compiled function that evaluates expr with args of the type specified.

compiledFunc = Compile[{x},
      Block[{sum = 1.0, inc = 1.0}, Do[inc = inc*x/i; sum = sum + inc, {i, 10000}]; sum]];

This is significantly faster than before,

Map[compiledFunc, arg]; // AbsoluteTiming

However, the performance can be improved yet further by allowing the function to list itself in a parallel manner (i.e. it no longer needs to be mapped):

compiledParallelFunc = Compile[{x},
      Block[{sum = 1.0, inc = 1.0}, Do[inc = inc*x/i; sum = sum + inc, {i, 10000}]; sum],
      RuntimeAttributes -> {Listable}, Parallelization -> True];

Any further optimisation would require additional effort:

compiledParallelFunc[arg]; // AbsoluteTiming

Note that a number of Wolfram Language functions automatically compile their arguments, or else cache results.

Symbolic C Generation

Should your code be succesfully compilable then you may also call your C Compiler to build a .DLL that can be run directly on your CPU rather than the Mathematica virtual machine:

compiledCCode = Compile[{x},
      Block[{sum = 1.0, inc = 1.0}, Do[inc = inc*x/i; sum = sum + inc, {i, 10000}]; sum],
      RuntimeAttributes -> {Listable}, Parallelization -> True, CompilationTarget -> "C"];

compiledCCode[arg]; // AbsoluteTiming

Internal`Bag` to replace Reap and Sow

Reap[Do[Sow[RandomInteger[{0, 10}]], {i, 100000}]]; // AbsoluteTiming

Reap and Sow do not work with Compile

cfFail = Compile[{}, Module[{},
        Reap[Do[Sow[RandomInteger[{0, 10}]], {i, 100000}]]]];

cfFail[]; // AbsoluteTiming

We can however use Bag, StuffBag and BagPart to the same effect.

cfBag = Compile[{}, Module[{b = Internal`Bag[{1}, 1]},
        Do[Internal`StuffBag[b, {RandomInteger[{0, 10}]}, 1], {i, 100000}];
        Internal`BagPart[b, All]]];

(x = cfz[]); // AbsoluteTiming

But there is Built in function for this

RandomInteger[{0, 10}, {100000}]; // AbsoluteTiming

Preventing Unnecessary Computations (1)

Optimising Plot

Plot has the HoldAll attribute, meaning that evaluations may be unnecessarily slow if an unevaluated expression is provided in the first argument:

Attributes[Plot]

Plot[
    Integrate[1/(1 + x)^2, {x, 0, n}],
    {n, 1, 15}] // AbsoluteTiming

If the evaluated expression is provided directly then Plot will be much faster, or else the function Evalaute could be wrapped around the first argument:

Plot[
    Evaluate[Integrate[1/(1 + x)^2, {x, 0, n}]],
    {n, 1, 15}] // AbsoluteTiming

There is another trick using the function With

With[
    {arg = Integrate[1/(1 + x)^2, {x, 0, n}]},
    Plot[arg, {n, 1, 15}]
    ] // AbsoluteTiming

Preventing Unnecessary Computations (2)

Dynamic content in Mathematica usually controls its updating in a sensible way, but there can be unexpected behaviour which you need to control to prevent unnecessary computations.

Consider this graphic:

ptdata = RandomReal[1, {10000, 2}];
Graphics[
  {Point /@ ptdata, PointSize[0.05], Red, Point[{0.5, 0.5}]}]

Suppose we wanted to build a Manipulate to move the red dot about, this would be the most obvious way to write it:

Manipulate[Graphics[
    {Point /@ ptdata, PointSize[0.05], Red, Point[{0.5, t}]}],
  {t, 0.05, 0.95}, SaveDefinitions -> True]

But it is very slow, because each time the variable t is changed the entire expression in the first argument is re-evaluated; the entire Graphics is re-drawn.

Many users don’t realise that the detailed language of DynamicModule can be utilised within Manipulate, the re-draw can be limited to just the red dot:

Manipulate[Graphics[
    {Point /@ ptdata, PointSize[0.05], Red, Dynamic[Point[{0.5, t}]]}],
  {t, 0.05, 0.95}, SaveDefinitions -> True]

Thank you for Listening!